In [1]:
#importing libraries
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns
from tqdm import tqdm_notebook
from igraph import Graph
import igraph
import Neuroscience_Project_modules as npm
import bct
import random
1. Connectivity graph¶
In [2]:
##paths
open_eyes = "./data/S004R01.edf"
closed_eyes = "./data/S004R02.edf"
In [3]:
freq_alpha = 10 ###10Hz Alpha_Wave_Frequency
1.1. (mandatory)¶
Estimate functional brain connectivity among 64 channels using one of the MVAR estimators: Partial Directed Coherence (PDC), Direct Transfer Function (DTF). Select one relevant frequency value. Apply a threshold so that the resulting binary connectivity matrices have network density equal to 20%. Create a graphical representation of the binary adjacency matrix.
Open Eyes
In [4]:
PDC_open = npm.Connectivity_Graph('PDC')
PDC_open.import_data(open_eyes)
PDC_open.connectivity(freq=freq_alpha,threshold=0.2,plot=True)
print(PDC_open.binary_adjacency_matrix)
plt.imshow(PDC_open.binary_adjacency_matrix,cmap='Greys',interpolation='nearest')
Out[4]:
Closed Eyes
In [5]:
PDC_closed = npm.Connectivity_Graph('PDC')
PDC_closed.import_data(closed_eyes)
PDC_closed.connectivity(freq=freq_alpha,threshold=0.2,plot=True)
print(PDC_closed.binary_adjacency_matrix)
plt.imshow(PDC_closed.binary_adjacency_matrix,cmap='Greys',interpolation='nearest')
Out[5]:
1.2. (class ‘A’)¶
Perform task 1.1 using both estimators (PDC and DTF).
Open Eyes
In [6]:
DTF_open = npm.Connectivity_Graph('DTF')
DTF_open.import_data(open_eyes)
DTF_open.connectivity(freq=freq_alpha,threshold=0.2,plot=True)
print(DTF_open.binary_adjacency_matrix)
plt.imshow(DTF_open.binary_adjacency_matrix,cmap='Greys',interpolation='nearest')
Out[6]:
Closed Eyes
In [7]:
DTF_closed = npm.Connectivity_Graph('DTF')
DTF_closed.import_data(closed_eyes)
DTF_closed.connectivity(freq=freq_alpha,threshold=0.2,plot=True)
print(DTF_closed.binary_adjacency_matrix)
plt.imshow(DTF_closed.binary_adjacency_matrix,cmap='Greys',interpolation='nearest')
Out[7]:
In [8]:
#Models Order computed from Akaike criterion
order_open = 7
order_closed = 6
1.3.(class ‘A’)¶
Perform task 1.1 using thresholds yielding the following density values: 1%, 5%, 10%, 20%, 30%, 50%.
In [9]:
m_open = []
m_closed = []
densities = [0.01,0.05,0.1,0.2,0.3,0.5]
for d in tqdm_notebook(densities):
PDC_open = npm.Connectivity_Graph('PDC')
PDC_open.import_data(open_eyes)
PDC_open.connectivity(freq=freq_alpha,threshold=d,order=order_open)
m_open.append(PDC_open.binary_adjacency_matrix)
PDC_closed = npm.Connectivity_Graph('PDC')
PDC_closed.import_data(closed_eyes)
PDC_closed.connectivity(freq=freq_alpha,threshold=d,order=order_closed)
m_closed.append(PDC_closed.binary_adjacency_matrix)
m = [m_open,m_closed]
In [10]:
cols = ['Case: {}'.format(col) for col in ['open','closed']]
rows = ['Density: {}'.format(row) for row in densities]
fig, axes = plt.subplots(nrows=6, ncols=2, figsize=(20,12))
for ax, col in zip(axes[0], cols):
ax.set_title(col)
for ax, row in zip(axes[:,0], rows):
ax.set_ylabel(row, rotation=90, size='large')
r = 0
for row in axes:
c = 0
for col in row:
col.imshow(m[c][r],cmap='Greys',interpolation='nearest')
c += 1
r += 1
fig.tight_layout()
fig.subplots_adjust(left=0.05, top=0.95)
plt.show()
1.4.(class ‘D’)¶
Considering the subset of 19 channels suggested in Figure 1 and Table 2, estimate the connectivity using PDC or DTF and apply a statistical validation method (asymptotic statistics 7 , resampling procedure 8 ,...) to filter out values that are not significantly different from 0 (PDC(i, j) ≠ 0 with p < 5%).
In [11]:
channels = 'Fp1 Fp2 F7 F3 Fz F4 F8 T7 C3 Cz C4 T8 P7 P3 Pz P4 P8 O1 O2'.split(' ')
Open Eyes
In [12]:
PDC_open = npm.Connectivity_Graph('PDC')
PDC_open.import_data(open_eyes,channels=channels)
PDC_open.connectivity(freq=freq_alpha,order=order_open)
PDC_open.significance(signf_threshold=0.05,channels=channels,visual=False,path=open_eyes,Nrep=1000)
#PDC_open.significance_matrix
In [13]:
plt.imshow(PDC_open.significance_matrix,cmap='jet_r')
Out[13]:
In [14]:
## associated binary adjacency matrix after insignificant connections removal
plt.imshow(PDC_open.binary_adjacency_matrix,cmap='Greys',interpolation='nearest')
Out[14]:
Closed Eyes
In [15]:
PDC_closed = npm.Connectivity_Graph('PDC')
PDC_closed.import_data(closed_eyes,channels=channels)
PDC_closed.connectivity(freq=freq_alpha,order=order_closed)
PDC_closed.significance(signf_threshold=0.05,channels=channels,visual=False,path=open_eyes,Nrep=1000)
#PDC_closed.significance_matrix
In [16]:
plt.imshow(PDC_closed.significance_matrix,cmap='jet_r')
Out[16]:
In [17]:
## associated binary adjacency matrix after insignificant connections removal
plt.imshow(PDC_closed.binary_adjacency_matrix,cmap='Greys',interpolation='nearest')
Out[17]:
1.5.(class ‘C’)¶
Make a topographical representation of the networks (see example in Figure 2). Cartesian coordinates of planar representation of EEG channels are available in Table 3 (see also the file channel_locations.txt).
(the choice of this task is advised in the case of 19-channel networks and/or density ≤ 5%).
Open Eyes
In [18]:
PDC_open = npm.Connectivity_Graph('PDC')
PDC_open.import_data(open_eyes)
PDC_open.connectivity(freq=freq_alpha,threshold=0.05,order=order_open)
PDC_open.show_graph('PDC_case_open')
Closed Eyes
In [19]:
PDC_closed = npm.Connectivity_Graph('PDC')
PDC_closed.import_data(closed_eyes)
PDC_closed.connectivity(freq=freq_alpha,threshold=0.05,order=order_closed)
PDC_closed.show_graph('PDC_case_closed')
1.6. (class ‘B’)¶
Perform task 1.1 considering a second frequency value belonging to a different EEG rhythm with respect to the first choice.
In [20]:
freq_beta = 20 ###20Hz Beta_Wave_Frequency
Open Eyes
In [21]:
PDC_open = npm.Connectivity_Graph('PDC')
PDC_open.import_data(open_eyes)
PDC_open.connectivity(freq=freq_beta,threshold=0.2,order=order_open)
print(PDC_open.binary_adjacency_matrix)
plt.imshow(PDC_open.binary_adjacency_matrix,cmap='Greys',interpolation='nearest')
Out[21]:
Closed Eyes
In [22]:
PDC_closed = npm.Connectivity_Graph('PDC')
PDC_closed.import_data(closed_eyes)
PDC_closed.connectivity(freq=freq_beta,threshold=0.2,order=order_closed)
print(PDC_closed.binary_adjacency_matrix)
plt.imshow(PDC_closed.binary_adjacency_matrix,cmap='Greys',interpolation='nearest')
Out[22]:
2. Graph theory indices¶
2.1. (mandatory)¶
Compute binary global (clustering coefficient, path length) and local (degree, in/out-degree) graph indices. List the highest 10 channels for local indices.
Open Eyes
In [23]:
IDX = npm.Graph_Theory_Indices("PDC")
IDX.import_data(open_eyes)
IDX.connectivity(freq=freq_alpha, threshold=0.2, plot=False,order=order_open)
print('\x1b[1;35m'+'\nGlobal indices'+'\x1b[0m \n')
print("APL (Average Path Length) and GCC (Global Clustering Coefficient)\n")
print('\x1b[1;39m'+"APL (w/o weights): "+'\x1b[0m', IDX.average_path_length(False))
print('\x1b[1;39m'+"GCC (w/o weights): "+'\x1b[0m', IDX.global_clustering_coefficient(False))
print("\n")
print('\x1b[1;35m'+'Local indices'+'\x1b[0m \n')
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m', IDX.local_indegree())
print('\x1b[1;39m'+"Local outdegree: "+'\x1b[0m', IDX.local_outdegree())
print('\x1b[1;39m'+"Local degree: "+'\x1b[0m', IDX.local_degree())
In [24]:
print('\x1b[1;35m'+'Top 10'+'\x1b[0m \n')
local_indegree_dict = IDX.local_indegree()
local_outdegree_dict = IDX.local_outdegree()
local_degree_dict = IDX.local_degree()
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m\n')
print(npm.top10(local_indegree_dict))
print('\n\x1b[1;39m'+"Local outdegree: "+'\x1b[0m\n')
print(npm.top10(local_outdegree_dict))
print('\n\x1b[1;39m'+"Local degree: "+'\x1b[0m\n')
print(npm.top10(local_degree_dict))
Closed Eyes
In [25]:
IDX = npm.Graph_Theory_Indices("PDC")
IDX.import_data(closed_eyes)
IDX.connectivity(freq=freq_alpha, threshold=0.2, plot=False,order=order_closed)
print('\x1b[1;35m'+'\nGlobal indices'+'\x1b[0m \n')
print("APL (Average Path Length) and GCC (Global Clustering Coefficient)\n")
print('\x1b[1;39m'+"APL (w/o weights): "+'\x1b[0m', IDX.average_path_length(False))
print('\x1b[1;39m'+"GCC (w/o weights): "+'\x1b[0m', IDX.global_clustering_coefficient(False))
print("\n")
print('\x1b[1;35m'+'Local indices'+'\x1b[0m \n')
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m', IDX.local_indegree())
print('\x1b[1;39m'+"Local outdegree: "+'\x1b[0m', IDX.local_indegree())
print('\x1b[1;39m'+"Local degree: "+'\x1b[0m', IDX.local_indegree())
In [26]:
print('\x1b[1;35m'+'Top 10'+'\x1b[0m \n')
local_indegree_dict = IDX.local_indegree()
local_outdegree_dict = IDX.local_outdegree()
local_degree_dict = IDX.local_degree()
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m\n')
print(npm.top10(local_indegree_dict))
print('\n\x1b[1;39m'+"Local outdegree: "+'\x1b[0m\n')
print(npm.top10(local_outdegree_dict))
print('\n\x1b[1;39m'+"Local degree: "+'\x1b[0m\n')
print(npm.top10(local_degree_dict))
2.2. (class ‘D’)¶
Search in the literature a definition of small-worldness index (i.e. an index describing the small-world organization of a network) and compute it.
Open Eyes
In [27]:
IDX = npm.Graph_Theory_Indices("PDC")
IDX.import_data(open_eyes)
IDX.connectivity(freq=freq_alpha, threshold=0.2, plot=False,order=order_open)
C = IDX.global_clustering_coefficient(False)
L = IDX.average_path_length(False)
print('\x1b[1;35m'+'\nAverage Path Length and Global Clustering Coefficient for\nGraph\n' + '\x1b[0m')
print('\x1b[1;39m'+ 'L '+'\x1b[0m' + '= ' + str(L),'\x1b[1;39m'+ ' C '+'\x1b[0m' + '= ' + str(C) )
APL_rnd_list = []
GCC_rnd_list = []
for cnt in range(10):
G_r_Adj = np.array([r for r in IDX.G.get_adjacency()])
G_r_Adj = bct.randmio_dir(G_r_Adj, 10)[0]
G_r_Adj = [list(r) for r in G_r_Adj]
G_r_Adj = [[G_r_Adj[i][j].item() for j in range(len(G_r_Adj))] for i in range(len(G_r_Adj))]
G_r = Graph.Adjacency(G_r_Adj)
APL_rnd_list.append(npm.average_path_length(G_r, False))
GCC_rnd_list.append(npm.global_clustering_coefficient(G_r, False))
L_r = np.mean(APL_rnd_list)
C_r = np.mean(GCC_rnd_list)
print('\x1b[1;35m'+'\nAverage Path Length and Global Clustering Coefficient for\nEquivalent Random Graph\n' + '\x1b[0m')
print('\x1b[1;39m'+ 'L_r '+'\x1b[0m' + '= ' + str(L_r),'\x1b[1;39m'+ ' C_r '+'\x1b[0m' + '= ' + str(C_r) )
APL_ltc_list = []
GCC_ltc_list = []
for cnt in range(10):
G_l_Adj = np.array([r for r in IDX.G.get_adjacency()])
G_l_Adj = bct.latmio_dir(G_l_Adj, 10)[0]
G_l_Adj = [list(r) for r in G_l_Adj]
G_l_Adj = [[G_l_Adj[i][j].item() for j in range(len(G_r_Adj))] for i in range(len(G_r_Adj))]
G_l = Graph.Adjacency(G_l_Adj)
APL_ltc_list.append(npm.average_path_length(G_l, False))
GCC_ltc_list.append(npm.global_clustering_coefficient(G_l, False))
L_l = np.mean(APL_ltc_list)
C_l = np.mean(GCC_ltc_list)
print('\x1b[1;35m'+'\nAverage Path Length and Global Clustering Coefficient for\nEquivalent Lattice Graph\n' + '\x1b[0m')
print('\x1b[1;39m'+ 'L_l '+'\x1b[0m' + '= ' + str(L_l),'\x1b[1;39m'+ ' C_l '+'\x1b[0m' + '= ' + str(C_l) )
SWI = (L - L_l)/(L_r - L_l) * (C - C_r)/(C_l - C_r)
print('\x1b[1;35m'+'\nSmall World Index\n' + '\x1b[0m')
print('\x1b[1;39m'+ 'SWI '+'\x1b[0m' + '=', SWI)
Closed Eyes
In [28]:
IDX = npm.Graph_Theory_Indices("PDC")
IDX.import_data(closed_eyes)
IDX.connectivity(freq=freq_alpha, threshold=0.2, plot=False,order=order_closed)
C = IDX.global_clustering_coefficient(False)
L = IDX.average_path_length(False)
print('\x1b[1;35m'+'\nAverage Path Length and Global Clustering Coefficient for\nGraph\n' + '\x1b[0m')
print('\x1b[1;39m'+ 'L '+'\x1b[0m' + '= ' + str(L),'\x1b[1;39m'+ ' C '+'\x1b[0m' + '= ' + str(C) )
APL_rnd_list = []
GCC_rnd_list = []
for cnt in range(10):
G_r_Adj = np.array([r for r in IDX.G.get_adjacency()])
G_r_Adj = bct.randmio_dir(G_r_Adj, 10)[0]
G_r_Adj = [list(r) for r in G_r_Adj]
G_r_Adj = [[G_r_Adj[i][j].item() for j in range(len(G_r_Adj))] for i in range(len(G_r_Adj))]
G_r = Graph.Adjacency(G_r_Adj)
APL_rnd_list.append(npm.average_path_length(G_r, False))
GCC_rnd_list.append(npm.global_clustering_coefficient(G_r, False))
L_r = np.mean(APL_rnd_list)
C_r = np.mean(GCC_rnd_list)
print('\x1b[1;35m'+'\nAverage Path Length and Global Clustering Coefficient for\nEquivalent Random Graph\n' + '\x1b[0m')
print('\x1b[1;39m'+ 'L_r '+'\x1b[0m' + '= ' + str(L_r),'\x1b[1;39m'+ ' C_r '+'\x1b[0m' + '= ' + str(C_r) )
APL_ltc_list = []
GCC_ltc_list = []
for cnt in range(10):
G_l_Adj = np.array([r for r in IDX.G.get_adjacency()])
G_l_Adj = bct.latmio_dir(G_l_Adj, 10)[0]
G_l_Adj = [list(r) for r in G_l_Adj]
G_l_Adj = [[G_l_Adj[i][j].item() for j in range(len(G_r_Adj))] for i in range(len(G_r_Adj))]
G_l = Graph.Adjacency(G_l_Adj)
APL_ltc_list.append(npm.average_path_length(G_l, False))
GCC_ltc_list.append(npm.global_clustering_coefficient(G_l, False))
L_l = np.mean(APL_ltc_list)
C_l = np.mean(GCC_ltc_list)
print('\x1b[1;35m'+'\nAverage Path Length and Global Clustering Coefficient for\nEquivalent Lattice Graph\n' + '\x1b[0m')
print('\x1b[1;39m'+ 'L_l '+'\x1b[0m' + '= ' + str(L_l),'\x1b[1;39m'+ ' C_l '+'\x1b[0m' + '= ' + str(C_l) )
SWI = (L - L_l)/(L_r - L_l) * (C - C_r)/(C_l - C_r)
print('\x1b[1;35m'+'\nSmall World Index\n' + '\x1b[0m')
print('\x1b[1;39m'+ 'SWI '+'\x1b[0m' + '=', SWI)
2.3. (class ‘B’)¶
Compare the global indices extracted from PDC and DTF connectivity estimations.
Open Eyes
In [29]:
IDX = npm.Graph_Theory_Indices("DTF")
IDX.import_data(open_eyes)
IDX.connectivity(freq=freq_alpha, threshold=0.2, plot=False,order=order_open)
print('\x1b[1;35m'+'\nGlobal indices'+'\x1b[0m \n')
print("APL (Average Path Length) and GCC (Global Clustering Coefficient)\n")
print('\x1b[1;39m'+"APL (w/o weights): "+'\x1b[0m', IDX.average_path_length(False))
print('\x1b[1;39m'+"GCC (w/o weights): "+'\x1b[0m', IDX.global_clustering_coefficient(False))
print("\n")
print('\x1b[1;35m'+'Local indices'+'\x1b[0m \n')
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m', IDX.local_indegree())
print('\x1b[1;39m'+"Local outdegree: "+'\x1b[0m', IDX.local_outdegree())
print('\x1b[1;39m'+"Local degree: "+'\x1b[0m', IDX.local_degree())
In [30]:
print('\x1b[1;35m'+'Top 10'+'\x1b[0m \n')
local_indegree_dict = IDX.local_indegree()
local_outdegree_dict = IDX.local_outdegree()
local_degree_dict = IDX.local_degree()
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m\n')
print(npm.top10(local_indegree_dict))
print('\n\x1b[1;39m'+"Local outdegree: "+'\x1b[0m\n')
print(npm.top10(local_outdegree_dict))
print('\n\x1b[1;39m'+"Local degree: "+'\x1b[0m\n')
print(npm.top10(local_degree_dict))
Closed Eyes
In [31]:
IDX = npm.Graph_Theory_Indices("DTF")
IDX.import_data(closed_eyes)
IDX.connectivity(freq=freq_alpha, threshold=0.2, plot=False,order=order_closed)
print('\x1b[1;35m'+'\nGlobal indices'+'\x1b[0m \n')
print("APL (Average Path Length) and GCC (Global Clustering Coefficient)\n")
print('\x1b[1;39m'+"APL (w/o weights): "+'\x1b[0m', IDX.average_path_length(False))
print('\x1b[1;39m'+"GCC (w/o weights): "+'\x1b[0m', IDX.global_clustering_coefficient(False))
print("\n")
print('\x1b[1;35m'+'Local indices'+'\x1b[0m \n')
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m', IDX.local_indegree())
print('\x1b[1;39m'+"Local outdegree: "+'\x1b[0m', IDX.local_outdegree())
print('\x1b[1;39m'+"Local degree: "+'\x1b[0m', IDX.local_degree())
In [32]:
print('\x1b[1;35m'+'Top 10'+'\x1b[0m \n')
local_indegree_dict = IDX.local_indegree()
local_outdegree_dict = IDX.local_outdegree()
local_degree_dict = IDX.local_degree()
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m\n')
print(npm.top10(local_indegree_dict))
print('\n\x1b[1;39m'+"Local outdegree: "+'\x1b[0m\n')
print(npm.top10(local_outdegree_dict))
print('\n\x1b[1;39m'+"Local degree: "+'\x1b[0m\n')
print(npm.top10(local_degree_dict))
2.4.(class ‘C’)¶
Study the behaviour of global graph indices in function of network density (see point 2.3 for density values).
(the choice of this task is advised in the case of selection of task 1.3).
Open Eyes
In [33]:
random.seed (7890)
APL_list = []
GCC_list = []
densities = [0.01,0.05,0.1,0.2,0.3,0.5]
for d in tqdm_notebook(densities):
IDX = npm.Graph_Theory_Indices("PDC")
IDX.import_data(open_eyes)
IDX.connectivity(freq=freq_alpha, threshold=d, plot=False,order=order_open)
APL_list.append(IDX.average_path_length(False))
GCC_list.append(IDX.global_clustering_coefficient(False))
In [34]:
plt.plot(densities, APL_list)
Out[34]:
In [35]:
plt.plot(densities, GCC_list)
Out[35]:
Closed Eyes
In [36]:
APL_list = []
GCC_list = []
densities = [0.01,0.05,0.1,0.2,0.3,0.5]
for d in tqdm_notebook(densities):
IDX = npm.Graph_Theory_Indices("PDC")
IDX.import_data(closed_eyes)
IDX.connectivity(freq=freq_alpha, threshold=d, plot=False,order=order_closed)
APL_list.append(IDX.average_path_length(False))
GCC_list.append(IDX.global_clustering_coefficient(False))
In [37]:
plt.plot(densities, APL_list)
Out[37]:
In [38]:
plt.plot(densities, GCC_list)
Out[38]:
2.5.(class ‘B’)¶
Make a topographical representation of local indices.
In [39]:
def plot_graph(G, name, degree_list):
G.es['arrow_size'] = [0.1 for edge in G.es]
visual_style = {}
visual_style["vertex_size"] = degree_list*10
visual_style["vertex_color"] = "red"
visual_style["vertex_label"] = G.vs["label"]
visual_style["edge_width"] = 0.2
visual_style["layout"] = G.vs["coords"]
graph = igraph.plot(G, bbox=(0, 0, 600, 600), **visual_style)
graph.save(name + '.png')
return(graph)
Open Eyes
In [40]:
IDX = npm.Graph_Theory_Indices("PDC")
IDX.import_data(open_eyes)
IDX.connectivity(freq=freq_alpha, threshold=0.2, plot=False,order=order_open)
In [41]:
degree_list = list(IDX.local_indegree().values())
plot_graph(IDX.G, "EO_indegree_graph", degree_list)
Out[41]:
In [42]:
degree_list = list(IDX.local_outdegree().values())
plot_graph(IDX.G, "EO_outdegree_graph", degree_list)
Out[42]:
In [43]:
degree_list = list(IDX.local_degree().values())
plot_graph(IDX.G, "EO_degree_graph", degree_list)
Out[43]:
Closed Eyes
In [44]:
IDX = npm.Graph_Theory_Indices("PDC")
IDX.import_data(closed_eyes)
IDX.connectivity(freq=freq_alpha, threshold=0.2, plot=False,order=order_closed)
In [45]:
degree_list = list(IDX.local_indegree().values())
plot_graph(IDX.G, "EC_indegree_graph", degree_list)
Out[45]:
In [46]:
degree_list = list(IDX.local_outdegree().values())
plot_graph(IDX.G, "EC_outdegree_graph", degree_list)
Out[46]:
In [47]:
degree_list = list(IDX.local_degree().values())
plot_graph(IDX.G, "EC_degree_graph", degree_list)
Out[47]:
2.6. (class ‘B’)¶
Compare the networks obtained with the analysis 1.6 in terms of graph indices.
(the choice of this task is advised only in the case of selection of task 1.6)
Open Eyes
In [48]:
IDX = npm.Graph_Theory_Indices("PDC")
IDX.import_data(open_eyes)
IDX.connectivity(freq=freq_beta, threshold=0.2, plot=False,order=order_open)
print('\x1b[1;35m'+'\nGlobal indices'+'\x1b[0m \n')
print("APL (Average Path Length) and GCC (Global Clustering Coefficient)\n")
print('\x1b[1;39m'+"APL (w/o weights): "+'\x1b[0m', IDX.average_path_length(False))
print('\x1b[1;39m'+"GCC (w/o weights): "+'\x1b[0m', IDX.global_clustering_coefficient(False))
print("\n")
print('\x1b[1;35m'+'Local indices'+'\x1b[0m \n')
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m', IDX.local_indegree())
print('\x1b[1;39m'+"Local outdegree: "+'\x1b[0m', IDX.local_outdegree())
print('\x1b[1;39m'+"Local degree: "+'\x1b[0m', IDX.local_degree())
In [49]:
print('\x1b[1;35m'+'Top 10'+'\x1b[0m \n')
local_indegree_dict = IDX.local_indegree()
local_outdegree_dict = IDX.local_outdegree()
local_degree_dict = IDX.local_degree()
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m\n')
print(npm.top10(local_indegree_dict))
print('\n\x1b[1;39m'+"Local outdegree: "+'\x1b[0m\n')
print(npm.top10(local_outdegree_dict))
print('\n\x1b[1;39m'+"Local degree: "+'\x1b[0m\n')
print(npm.top10(local_degree_dict))
Closed Eyes
In [50]:
IDX = npm.Graph_Theory_Indices("PDC")
IDX.import_data(closed_eyes)
IDX.connectivity(freq=freq_beta, threshold=0.2, plot=False,order=order_closed)
print('\x1b[1;35m'+'\nGlobal indices'+'\x1b[0m \n')
print("APL (Average Path Length) and GCC (Global Clustering Coefficient)\n")
print('\x1b[1;39m'+"APL (w/o weights): "+'\x1b[0m', IDX.average_path_length(False))
print('\x1b[1;39m'+"GCC (w/o weights): "+'\x1b[0m', IDX.global_clustering_coefficient(False))
print("\n")
print('\x1b[1;35m'+'Local indices'+'\x1b[0m \n')
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m', IDX.local_indegree())
print('\x1b[1;39m'+"Local outdegree: "+'\x1b[0m', IDX.local_outdegree())
print('\x1b[1;39m'+"Local degree: "+'\x1b[0m', IDX.local_degree())
In [51]:
print('\x1b[1;35m'+'Top 10'+'\x1b[0m \n')
local_indegree_dict = IDX.local_indegree()
local_outdegree_dict = IDX.local_outdegree()
local_degree_dict = IDX.local_degree()
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m\n')
print(npm.top10(local_indegree_dict))
print('\n\x1b[1;39m'+"Local outdegree: "+'\x1b[0m\n')
print(npm.top10(local_outdegree_dict))
print('\n\x1b[1;39m'+"Local degree: "+'\x1b[0m\n')
print(npm.top10(local_degree_dict))
2.7. (class ‘C’)¶
Perform point 2.1 considering the weighted version of the graph indices definitions.
Open Eyes
In [52]:
IDX = npm.Graph_Theory_Indices("PDC")
IDX.import_data(open_eyes)
IDX.connectivity(freq=freq_alpha, threshold=0.2, plot=False,order=order_open)
print('\x1b[1;35m'+'\nGlobal indices'+'\x1b[0m \n')
print("APL (Average Path Length) and GCC (Global Clustering Coefficient)\n")
print('\x1b[1;39m'+"APL (w/ weights): "+'\x1b[0m', IDX.average_path_length(True))
print('\x1b[1;39m'+"GCC (w/ weights): "+'\x1b[0m', IDX.global_clustering_coefficient(True))
print("\n")
print('\x1b[1;35m'+'Local indices'+'\x1b[0m \n')
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m', IDX.local_indegree(True))
print('\x1b[1;39m'+"Local outdegree: "+'\x1b[0m', IDX.local_indegree(True))
print('\x1b[1;39m'+"Local degree: "+'\x1b[0m', IDX.local_indegree(True))
In [53]:
print('\x1b[1;35m'+'Top 10'+'\x1b[0m \n')
local_indegree_dict = IDX.local_indegree(True)
local_outdegree_dict = IDX.local_outdegree(True)
local_degree_dict = IDX.local_degree(True)
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m\n')
print(npm.top10(local_indegree_dict))
print('\n\x1b[1;39m'+"Local outdegree: "+'\x1b[0m\n')
print(npm.top10(local_outdegree_dict))
print('\n\x1b[1;39m'+"Local degree: "+'\x1b[0m\n')
print(npm.top10(local_degree_dict))
Closed Eyes
In [54]:
IDX = npm.Graph_Theory_Indices("PDC")
IDX.import_data(closed_eyes)
IDX.connectivity(freq=freq_alpha, threshold=0.2, plot=False,order=order_open)
print('\x1b[1;35m'+'\nGlobal indices'+'\x1b[0m \n')
print("APL (Average Path Length) and GCC (Global Clustering Coefficient)\n")
print('\x1b[1;39m'+"APL (w/ weights): "+'\x1b[0m', IDX.average_path_length(True))
print('\x1b[1;39m'+"GCC (w/ weights): "+'\x1b[0m', IDX.global_clustering_coefficient(True))
print("\n")
print('\x1b[1;35m'+'Local indices'+'\x1b[0m \n')
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m', IDX.local_indegree(True))
print('\x1b[1;39m'+"Local outdegree: "+'\x1b[0m', IDX.local_indegree(True))
print('\x1b[1;39m'+"Local degree: "+'\x1b[0m', IDX.local_indegree(True))
In [55]:
print('\x1b[1;35m'+'Top 10'+'\x1b[0m \n')
local_indegree_dict = IDX.local_indegree(True)
local_outdegree_dict = IDX.local_outdegree(True)
local_degree_dict = IDX.local_degree(True)
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m\n')
print(npm.top10(local_indegree_dict))
print('\n\x1b[1;39m'+"Local outdegree: "+'\x1b[0m\n')
print(npm.top10(local_outdegree_dict))
print('\n\x1b[1;39m'+"Local degree: "+'\x1b[0m\n')
print(npm.top10(local_degree_dict))
3. Motif analysis¶
3.1 (mandatory)¶
Perform motifs analysis to investigate the presence of 3-node configurations in the networks: determine their frequency and statistical significance (motifs, anti-motifs).
In [56]:
## alpha waves
## threshold: 0.05
## method: PDC
## method for p: AIC
freq=10
threshold=0.05
method = 'PDC'
Open Eyes
In [57]:
# Create EO connectivity graph for 10Hz frequence
# then save the edges in txt file for mfinder app
EO = npm.Graph_for_motifs(method, 'EO', density_threshold=threshold, freq=freq, folder='data/')
EO.connectivity(freq=freq,threshold=threshold,plot=True)
EO.preprocess_for_mfinder()
In [58]:
# performing motifs analysis with mfinder:
# mfinder ../data/EO_10Hz_0.05dt_Motifs.txt -s 3 -r 2000 -nu
# mfinder ../data/EC_10Hz_0.05dt_Motifs.txt -s 3 -r 2000 -nu
In [59]:
with open('data/EO_10Hz_0.05dt_Motifs_OUT.txt', 'r') as f:
print (f.read())
Summary for EO:
===============
2 motifs detected:
| ID | 6, 108 |
===============
1 anti-motifs detected:
| ID | x |
Closed Eyes
In [60]:
# Create EC connectivity graph for 10Hz frequence
# then save the edges in txt file for mfinder app
EC = npm.Graph_for_motifs(method, 'EC', density_threshold=threshold, freq=freq_alpha, folder='data/')
EC.connectivity(freq=freq,threshold=threshold,plot=True)
EC.preprocess_for_mfinder()
In [61]:
with open('data/EC_10Hz_0.05dt_Motifs_OUT.txt', 'r') as f:
print (f.read())
Summary for EC:
===============
4 motifs detected:
| ID | 38, 46, 108, 238|
===============
x anti-motifs detected
| ID | x |
3.2. (class ‘C’)¶
For the motif with pattern A → B ← C, create a topographical representation of the networks considering only the connections involved in this configuration.
Open Eyes
In [62]:
# Motif file structure: EO
with open('data/EO_10Hz_0.05dt_Motifs_MEMBERS.txt', 'r') as f:
for _ in range(9):
print(f.readline())
print('\t...')
In [63]:
EO.plot_motifs_36('data/EO_10Hz_0.05dt_Motifs_MEMBERS.txt',
name='data/EO_0.05_motif36', save=False)
Out[63]:
Closed Eyes
In [64]:
# Motif file structure: EC
with open('data/EC_10Hz_0.05dt_Motifs_MEMBERS.txt', 'r') as f:
for _ in range(9):
print(f.readline())
print('\t...')
In [65]:
EC.plot_motifs_36('data/EC_10Hz_0.05dt_Motifs_MEMBERS.txt',
name='data/EC_0.05_motif36', save=False)
Out[65]:
3.3. (class ‘C’)¶
Choose a channel selected in parieto-occipital scalp region and determine the motifs which involve it.
Open Eyes
In [66]:
df_prova = EO.find_motifs(channel='Po3')
df_prova
Out[66]:
Closed Eyes
In [67]:
EC.find_motifs('Po3')
Out[67]:
Additional: let's see the motifs in configuration A-->B<--C involving Po3 channel¶
In [68]:
# channel selected: Po3 (57)
# in our notation (from 0 to 64) Po3 correspond to 56
# we can directly plot it with our figure
# remaking member analysis
In [69]:
# plot for EO (pattern A->B<-C)
EO.plot_motifs_36('data/EO_10Hz_0.05dt_Motifs_MEMBERS.txt', 'EO_motifsPlot_Po3',
single_channel=56, motifcolor="red", save=False, tipsize=0.5)
Out[69]:
In [70]:
# plot for EC (pattern A->B<-C)
EC.plot_motifs_36('data/EC_10Hz_0.05dt_Motifs_MEMBERS.txt', 'EC_motifsPlot_Po3',
single_channel=56, motifcolor="red", save=False, tipsize=0.5)
Out[70]:
4. (class ‘E’)¶
Perform the same analysis described in task 3.1 considering 4-node motifs
The same analysis has been performed using mfinder. The netork used was the same of point 3.1
Open Eyes
In [71]:
with open('data/EO_10Hz_0.05dt_4_Motifs_OUT.txt', 'r') as f:
for _ in range(51):
print(f.readline())
print('...\n...\nCheck EO_10Hz_0.05dt_4_Motifs_OUT.txt for more info')
Closed Eyes
In [72]:
with open('data/EC_10Hz_0.05dt_4_Motifs_OUT.txt', 'r') as f:
for _ in range(51):
print(f.readline())
print('...\n...\nCheck EC_10Hz_0.05dt_4_Motifs_OUT.txt for more info')
4. Community detection¶
4.1. (mandatory)¶
Determine number and composition (i.e. list of nodes) of the communities obtained applying one of the algorithms introduced during the course.
Louvain_Algorithm - Optimal_Modularity¶
Open Eyes
In [73]:
PDC_open = npm.Connectivity_Graph('PDC')
PDC_open.import_data(open_eyes)
PDC_open.connectivity(freq=freq_alpha,threshold=0.2,plot=True)
Closed Eyes
In [74]:
PDC_closed = npm.Connectivity_Graph('PDC')
PDC_closed.import_data(closed_eyes)
PDC_closed.connectivity(freq=freq_alpha,threshold=0.2,plot=True)
4.2. (class ‘B’)¶
Make a graphical representation of the community structure in both rest conditions.
Open Eyes
In [75]:
comms_mod_open,graph = PDC_open.optimal_modularity_community_detection(name='optimal_modularity_open_eyes')
graph
Out[75]:
Closed Eyes
In [76]:
comms_mod_closed,graph = PDC_closed.optimal_modularity_community_detection(name='optimal_modularity_closed_eyes')
graph
Out[76]:
4.3. (class ‘C’)¶
Compare the community structure obtained by means of two different methods (modularity-based vs information theory-based approaches).
Infomap_Algorithm_random_walk - Information_Theory¶
Open Eyes
In [77]:
comms_info_open,graph = PDC_open.infomap_community_detection(name='information_theory_open_eyes')
graph
Out[77]:
Communities comparison for open_eyes_case¶
In [78]:
M = np.zeros([len(comms_mod_open),len(comms_info_open)])
for i in range(len(comms_mod_open)):
for j in range(len(comms_info_open)):
#JACCARD SIMILARITY
M[i,j] = len(set(comms_mod_open[i]).intersection(set(comms_info_open[j])))/len(set(comms_mod_open[i]).union(set(comms_info_open[j])))
M = pd.DataFrame(M)
plt.figure(figsize=(9,9))
ax = sns.heatmap(M,xticklabels=M.columns+1,yticklabels=M.columns+1,annot=True,cmap='viridis',linewidths=1,square=True)
ax.invert_yaxis()
plt.xlabel('Infomap',fontsize=12)
xlabels = comms_info_open.keys()
ylabels = comms_mod_open.keys()
plt.xticks(np.arange(len(comms_info_open))+0.5, xlabels,)
plt.yticks(np.arange(len(comms_mod_open))+0.5, ylabels)
plt.ylim(0,len(comms_mod_open))
plt.ylabel('Optimal_Modularity_Method',fontsize=12)
plt.title('Jaccard Similarity',fontsize=15)
plt.tight_layout()
plt.show()
In [79]:
comms_info_closed,graph = PDC_closed.infomap_community_detection(name='information_theory_closed_eyes')
graph
Out[79]:
Communities comparison for closed_eyes_case¶
In [80]:
M = np.zeros([len(comms_mod_closed),len(comms_info_closed)])
for i in range(len(comms_mod_closed)):
for j in range(len(comms_info_closed)):
#JACCARD SIMILARITY
M[i,j] = len(set(comms_mod_closed[i]).intersection(set(comms_info_closed[j])))/len(set(comms_mod_closed[i]).union(set(comms_info_closed[j])))
M = pd.DataFrame(M)
plt.figure(figsize=(9,9))
ax = sns.heatmap(M,xticklabels=M.columns+1,yticklabels=M.columns+1,annot=True,cmap='viridis',linewidths=1,square=True)
ax.invert_yaxis()
plt.xlabel('Infomap',fontsize=12)
xlabels = comms_info_closed.keys()
ylabels = comms_mod_closed.keys()
plt.xticks(np.arange(len(comms_info_closed))+0.5, xlabels,)
plt.yticks(np.arange(len(comms_mod_closed))+0.5, ylabels)
plt.ylim(0,len(comms_mod_closed))
plt.ylabel('Optimal_Modularity_Method',fontsize=12)
plt.title('Jaccard Similarity',fontsize=15)
plt.tight_layout()
plt.show()
plot of how number of communities varies for both resting cases for both methods¶
In [81]:
##densities threshold
densities = [0.05,0.1,0.2,0.3,0.5]
m1 = []
m2 = []
m3 = []
m4 = []
for d in tqdm_notebook(densities):
PDC_open = npm.Connectivity_Graph('PDC')
PDC_open.import_data(open_eyes)
PDC_open.connectivity(freq=freq_alpha,threshold=d,plot=False,order=order_open)
PDC_closed = npm.Connectivity_Graph('PDC')
PDC_closed.import_data(closed_eyes)
PDC_closed.connectivity(freq=freq_alpha,threshold=d,plot=False,order=order_closed)
comms_open_mod = PDC_open.optimal_modularity_community_detection(visual = False)
comms_closed_mod = PDC_closed.optimal_modularity_community_detection(visual = False)
comms_open_info = PDC_open.infomap_community_detection(visual = False)
comms_closed_info = PDC_closed.infomap_community_detection(visual = False)
m1.append(len(comms_open_mod))
m2.append(len(comms_open_info))
m3.append(len(comms_closed_mod))
m4.append(len(comms_closed_info))
Open Eyes
In [82]:
# set width of bar
barWidth = 0.25
# set height of bar
bars1 = m1
bars2 = m2
# Set position of bar on X axis
r1 = np.arange(len(bars1))
r2 = [x + barWidth for x in r1]
# Make the plot
plt.bar(r1, bars1, color='#7f6d5f', width=barWidth, edgecolor='white', label='Louvain')
plt.bar(r2, bars2, color='#557f2d', width=barWidth, edgecolor='white', label='Infomap')
# Add xticks on the middle of the group bars
plt.xlabel('graph density', fontweight='bold')
plt.xticks([r + barWidth for r in range(len(bars1))], densities)
# Create legend & Show graphic
plt.title('Number of communities - Open Eyes Case')
plt.legend()
plt.show()
Closed Eyes
In [83]:
# set width of bar
barWidth = 0.25
# set height of bar
bars1 = m3
bars2 = m4
# Set position of bar on X axis
r1 = np.arange(len(bars1))
r2 = [x + barWidth for x in r1]
# Make the plot
plt.bar(r1, bars1, color='#7f6d5f', width=barWidth, edgecolor='white', label='Louvain')
plt.bar(r2, bars2, color='#557f2d', width=barWidth, edgecolor='white', label='Infomap')
# Add xticks on the middle of the group bars
plt.xlabel('graph density', fontweight='bold')
plt.xticks([r + barWidth for r in range(len(bars1))], densities)
# Create legend & Show graphic
plt.title('Number of communities - Closed Eyes Case')
plt.legend()
plt.show()
